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Abstract 

The observed slow running of the gauge coupling in SU(3) lattice gauge theory with two flavors 
of color sextet fermions naturally suggests it is a theory with one relevant coupling, the fermion 
mass, and that at zero mass correlation functions decay algebraically. I perform a finite-size scaling 
study on simulation data at two values of the bare gauge coupling with this assumption and observe 
a common exponent for the scaling of the correlation length with the fermion mass, y m ~ 1.5. An 
analysis of the scaling of valence Dirac eigenvalues at one of these bare couplings produces a similar 
number. 
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I. INTRODUCTION AND THEORETICAL BACKGROUND 



Recently, many researchers |, |, |, 0, 0, 0, fi S 0, El El, H 0, El, EE E3, H E3 have 
begun to use lattice methods to study field theories which might be candidates for strongly 
coupled beyond - Standard Model new physics (23|. These models typically involve gauge 
fields and a large number of fermion degrees of freedom, either many flavors of fundamental 
representation fermions (where the discussion goes back to Refs. jzj. [25| ) or a smaller number 
of flavors of higher- dimensional representation fermions (strongly emphasized by Refs. [26|, 



23,1231). 

The usual description of renormalization for a gauge theory coupled to massless fermions 
classifies the possibilities for its behavior according to how the gauge coupling flows under 
rescaling to the infrared: it could flow to strong coupling in the case of a confining theory, 
or to zero, for a trivial theory, or to a fixed point at some nonzero value, g 2 *. The latter 
case is referred to as an infrared-attractive fixed point (IRFP) theory. This is a phase with 
no confinement, no chiral symmetry breaking, and algebraic decay of correlation functions. 
Evidence has been presented that several such theories exist: SU(3) gauge theory with 
Nf = 12 flavors of fundamental fermions @,J3, SU(3) gauge theory with Nf = 2 flavors of 
sextet fermions (the subject of this paper) [5[] and SU(2) gauge theory with Nf — 2 flavors 
of adjoint fermions (l5|. 

While for confining theories the gauge coupling is relevant (the Gaussian fixed point g = 
is unstable), that is not the case for IRFP theories. The distance of the bare gauge coupling 
from its fixed point value {g 2 —g 2 *) is an irrelevant coupling. The critical surface encompasses 
a wide range of values of bare gauge couplings. The fermion mass is a relevant coupling of 
an IRFP theory, since it must be fine-tuned at the UV scale to reach the critical surface (i.e. 
m q = 0). The situation is completely equivalent to that of an order-disorder transition in a 
magnetic system. The only difference is that the relevant direction is parameterized by the 
quark mass m q , rather than by the reduced temperature t = (T — T c )/T c of the magnet. (A 
closer analogy is to a system which has been fine tuned to its Curie point and then placed 
in an external magnetic field. The external field breaks the underlying global symmetry just 
as a quark mass explicitly breaks chiral symmetry.) [2J| 

The framework to describe the physics of these systems is standard. It involves a set of 
scaling operators which evolve independently and multiplicatively and whose renormalization 
group (RG) equations have linear zeroes; the RG equation for the change in the ith coupling, 
under a scale change by a factor s, linearized around its fixed-point value g*, is 

s ^ = Vi{9i- 9*)- (!) 

The mass is the relevant coupling and its exponent will be labeled y m . Tuning the mass to 
zero causes the correlation length to diverge algebraically, 

i_ 

f ~ m q ym . (2) 

The exponent y m is related to the evolution of the condensate in beyond-Standard Model 
phenomenology. Under a change of scale from s\ to S2, the condensate runs as 

f 32 dfx 

(i>i>)s 2 = (iH>)si ex P / — d i^) (3) 
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where d is its scaling dimension, (d = 3 + where 7^ is its anomalous dimension). 
Because the combination mxpip is an RG invariant, y m = 4 — d. 

Thus y m is an interesting object [3(j. The subject of this paper is a calculation of y m for 
SU(3) gauge fields coupled to Nf — 2 flavors of fermions in the sextet (color symmetric) 
representation. 

Most recent lattice work on candidate theories is devoted to answering the question of 
whether or not the gauge coupling runs to a fixed point. The analysis in this paper simply 
assumes that the gauge coupling runs very slowly. What happens in that case can be seen 
from an elaboration of Eq. [3J Suppose that we have a correlation function measured on 
a momentum scale which is large compared to other possible scales in the theory. Then, 
textbooks tell us that a generic correlation function (with engineering dimension d n and 
associated anomalous dimension 7) scales as 



When the integral is dominated by scales where the coupling is given by its fixed point value, 
then T(k) ~ k dn+1 where 7 = j(g*). The correlator has power law behavior. (This behavior 
is modified at short distances by non-scaling terms in the action.) 

If we have a real fixed point, then the critical exponents (in this case, y m ) will not depend 
on the value of the bare coupling. However, imagine that we do not have a fixed point. Then 
the integral in Eq. H] will depend on s. But if the change in the coupling over the range of 
the integral is small, 7(#(i)) will remain unchanged, and one will again observe a power law 
behavior for correlators. Of course, the value of the exponent 7 would change as the bare 
coupling were varied. Whether one is seeing a real exponent (a constant y m ) or an effective 
one (that is, one is mapping out y m (g)) can be tested by determining y m for several values 
of the bare gauge coupling. 

This discussion would be redundant if we knew from other sources whether the theory 
was in the IRFP phase or not. It is necessary to make it because this part of the story 
is not complete. It is quite clear that, in my system's weak coupling phase, observables 
(correlation lengths) depend quite strongly on the quark mass and only weakly on the gauge 
coupling. The correlation length saturates at a value proportional to the size of the system 
as the quark mass is tuned to zero. 

And in the weak coupling phase, the running coupling does run slowly. Last year S vet it- 
sky, Shamir and I [B| presented evidence that Nf = 2 sextet fermions and SU (3) gauge fields 
had an IRFP. Those simulations were done with a different bare action than the one used 
here. We are repeating and extending our calculation of the running coupling using the 
present bare action. That analysis is at present incomplete (3l| . At this point in time, we 
know that throughout the weak coupling phase the running coupling (defined through the 
Schrodinger Functional (SF) scheme) runs more slowly than two-loop perturbation theory 
predicts. This itself is slow running compared to the familiar case of Nf = 2 fundamen- 
tal representation QCD. That this should be so is obvious from the beta function, but the 
numerical values are worth mentioning. At the two bare couplings (3 = Q/g 2 where I will 
claim a measurement of y m , the SF coupling (measured on 6 4 volumes) is about g$ F = 2.5 at 
(3 = 5.2 and about 3.4 at (3 = 4.8. (Parameter sets will be given below.) The change of g^ F in 
a scale factor of 2 from integrating the two-loop perturbative beta function is about -0.2 and 
-0.3 respectively - about ten per cent. This justifies treating the gauge coupling as if it were 
not running and looking for deviations. For the more familiar case of SU(3) gauge group 




(4) 
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and Nf = 2 fundamental fermions, the change would be -0.8 and -2.0 at these couplings. 
(And in this case the observed coupling runs faster than the perturbative result. [32j].) 

Looking ahead to my answer, I claim that y m ~ 1.5 at the two values of the bare coupling 
where I measured it. 

I will investigate three different ways to measure y m .The successful ones employ finite 
size scaling arguments for the response of observables to the simulation volume. 

The first way uses the correlation length directly. When it grows to be the size of the 
simulation volume, Eq. [2] breaks down. Finite size scaling arguments allow us to use the 
breakdown to determine the exponent. This test is done at two values of the bare coupling. 

The two other tests involve working with extensive quantities. As m q is tuned to zero, 
the singular part of the free energy scales as 

f s (m q ) = mf^(A 1 + A 2 m)^ ym + ...). (5) 

where D is the system's dimensionality (here D = 4), and the prefactor is just (the 
correlation length provides the appropriate dimensional factor). A\ and A 2 are non- universal 
constants and yi is the biggest non-leading exponent. This is most likely the exponent y g of 
the gauge coupling, g 2 — g 2 *, which can be determined from the beta function as measured 
in (for example) Schrodinger functional simulations at m q = 0. It is probably small; in 
Ref. [29j], Hasenfratz and I estimated that it was close to zero in this and related theories. 
So the condensate, (ipip) = S(m g ), scales with m q as 



df s 



~ < (6) 



ql dm q q 

where a = D/y m — 1 . This is exactly like the relation of the specific heat exponent to the 
correlation length exponent at a conventional paramagnetic- ferromagnetic critical point. 

Working with the condensate has the problem that the particular lattice fermions used in 
this simulation, Wilson-type fermions, explicitly break chiral symmetry in the action. This 
introduces an additive shift to the condensate. So I will use partial quenching: I will take 
configurations generated with lattice fermions which do not have exact chiral symmetry, and 
measure the Dirac spectrum using valence quarks which are an implementation of lattice 
fermions with exact chiral symmetry (overlap fermions). One can question whether the 
valence fermion sees a faithful realization of what is happening in the equilibrium distribution 
of real dynamical variables. In usual (low- A/ fundamental QCD) this is believed to be the 
case. I think that for a first study, what I am going to do is adequate. (It would of course 
be better to do simulations with chiral dynamical fermions, as was done by the authors of 
Ref. but they are presently too expensive for simulations at large volume.) 

The condensate has UV-sensitive pieces j33[, (ijjip)uv ~ C\m q + C^m? q + . . . where C\ ~ 
1/a 2 and C3 ~ log a. This masks the non-analytic behavior. This is quite similar to the 



situation in finite temperature QCD, precisely at T = T c , where [34 



S(a, m q , T) ~ C\m q l a 2 + c$m^ s + analytic. (7) 

An IRFP theory is different from QCD in that there are no Goldstone bosons (which con- 
tribute their own non-analytic piece to the condensate, below T c ). It is also different from 
QCD in that while in QCD, Eq. [7J applies only at T c , in an IRFP theory Eq. [7] gives the 
behavior of the condensate throughout the basin of attraction of the IRFP. Unfortunately, 
for my data set, the UV terms dominate E(m). 
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A better way to attack y m through the condensate involves the Banks-Casher relational: 
between the condensate and the density of eigenvalues A of the Dirac operator p(A). At 
nonzero mass it is 

E( m ,, = -/p(AMA^. (8) 

If the massless theory is conformal, and the condensate S(m g ) ~ m° for small mass, then 
p(A) also scales as A a . A finite-size scaling argument 36[ relates the scaling for the density p 



to the scaling of the value of individual eigenvalues. If we consider the average value of the 
ith eigenvalue of the Dirac operator in a box of volume V = L D , and if p(A) ~ A a , then we 
expect 

(A,) ~ (tY (9) 



where the exponent is 



L 



"=TT-a (10) 



For the case of an IRFP theory, p is equal to y m , the leading exponent. 

(To derive this, note p ~ A" means that eigenvalues are uniformly distributed in an 
N = a + 1 dimensional space of volume V = R N , 



N 

IT 



A = ^E^) 1/2 n t = l^...R (11) 



i=l 



so an eigenvalue scales as Aj ~ 1/R = (l/V)* = (l/V)^ 1 . Now suppose we are in D 
physical dimensions; in a box of volume V, there are V = L D modes, from which Eq. [TU] is 
obtained.) 

One example of this formula is free field theory: a = D — 1 and p = 1. Another is the case 
of chiral symmetry breaking encoded in the usual formulas of its Random Matrix Theory 
analog: a = so p(A) — > p a constant, and p = D. This is (Aj) ~ 1/V, which is equivalent 
to the usual statement that the eigenvalue spectrum depends on the dimensionless product 
XEV. 

For the case of a system which exhibits chiral symmetry breaking, there is a tight theoret- 
ical description of the behavior of the lowest eigenvalues of the Dirac operator, which allows 
one to relate delicate features of the spectrum to the low energy constants of the theory (the 
condensate, the pseudoscalar decay constant, and possibly others). This description is based 
on Random Matrix Theory (RMT). However, if a system is in a chirally- restored phase, there 
are no longer RMT predictions to compare results against. Previous work shows that our 
target theory has a weak-coupling phase which is chirally-restored. The restoration of chiral 
symmetry is observed through regularities in the spectrum of screening masses as well as 
the behavior of the pseudoscalar decay constant as a function of quark mass. Therefore, 
the analysis reported here only uses the simplest property of the eigenvalues, namely their 
scaling with system size. 

The paper proceeds as follows. In Sec. [IT] I describe details of the simulations. In Sec. IIIIII 
compute the correlation length exponent using finite size scaling. Next, in Sec. [IV] I examine 
the mass dependence of the chiral condensate, and finally in Sec. |V]I perform a scaling test 
for eigenvalues of the valence Dirac operator. I conclude with a discussion of my results. 



In an earlier preprint 38[ I tried to do several of the analyses I report on, in this paper. 



In Ref. [Ill] we observed that in the deconfined phase of sextet QCD, the Polyakov loop 
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ordered in one of the negative real directions, roughly along one of the complex elements of 
Z(3), (Re(TrP(x)) < 0, Im(TrP(x)) 7^ 0). I believed that was the general situation for 
this theory and all the simulations were done in those vacua. After that paper appeared, 
detailed studies of the phase structure (at smaller volumes) by Machtey and Svetitsky [3^ 
showed that the true vacuum of the deconfmed phase is in fact the one in which the Polyakov 
loop is real and positive. All the results of Ref. [38[ only apply to metastable vacua. Their 
conclusion renders it too uninteresting to publish. 



II. NUMERICAL TECHNIQUES AND BACKGROUND 

I performed simulations on a system with SU(3) gauge fields and two flavors of dynamical 
fermions in the symmetric (sextet) representation of the color gauge group. The lattice 
action is defined by the single-plaquette gauge action and a Wilson fermion action with 
added clover term 3j|. The fermion action employs the differentiable hypercubic smeared 



link of Ref. 40j, from which the symmetric- representation gauge connection for the fermion 



operator is constructed. No tadpole-improvement is used and the clover coefficient is set 



to its tree- level value. The smearing parameters for the links are the same as in Ref. [40 
q;i = 0.75, ct2 = 0.6, 03 = 0.3. The bare parameters which are inputs to the simulation are 
the gauge coupling (3 = Q/g 2 and the fermion hopping parameter k. The integration is done 



with one additional heavy pseudo- fermion field as suggested by Hasenbusch [411 ]. multiple 
time scales [42|], and a second-order Omelyan integrator [43| . 

The routines for simulating sextet-representation fermions were developed with (and 
mostly by) B. Svetitsky and Y. Shamir. The dynamical fermion algorithm was adapted 
from a program written by A. Hasenfratz, R. Hoffmann and S. Schaefer [44| All computer 
code is based on the publicly available package of the MILC collaboration [451 ]. 

Simulation volumes range up to 16 4 sites, and typical data sets range from a few hundred 
to a thousand trajectories. I recorded lattices every five trajectories (of unit simulation time 
trajectory length) and collected 40 lattices per parameter set for the calculation of spectral 
observables and overlap eigenvalues. 

Correlation lengths are taken to be inverses of screening masses, from correlators of 
operators at different separations measured along one of the spatial directions of the lattice. 
The trick of combining periodic and anti-periodic boundary conditions for valence quarks 
(H, 47, H, 4^1 is used in these measurements. 



Throughout this work, instead of quoting k, I will use the the Axial Ward Identity (AWI) 
quark mass, defined through 

X X 

where Aq = ^7075^5 P — V^sVs an d X is any source. The derivative is taken to be the 
naive difference operator {d lx f{x) = {f{x + (xa) — f{x — jla))/ (2a)). For X I used a Coulomb 
gauge-fixed Gaussian source. 

Machtey and Svetitsky (3?J have performed careful studies of the phase structure of this 
theory and have shown that that the true vacuum (in the deconfmed phase) is the one 
in which the expectation value of the Polyakov loop is real and positive. All simulations 
reported here were done in this vacuum. 

The valence Dirac operator whose eigenvalues are used in Sec. |V] is the overlap operator 



50j, l5JJ. Details of the particular implementation of the action are described in Refs. [5J 
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53l . l54j . l55l . [56|. The only new ingredient is the application to symmetric-representation 
fermions, using the same combination of hypercubic link and projection into the fermionic 
representation as for the dynamical fermions. Eigenvalues of the squared Hermitian Dirac 
operator D^D are computed using the "Primme" package of McCombs and Stathopoulos[57| 
and split apart in the usual way. Eigenvalues are quoted after stereographic projection. 

There is a potential problem with this analysis involving the index theorem, relating the 
winding number of the gauge field Q to the number of Dirac fermion zero modes, 

index = 2T(R)Q. (13) 

T(R) is the Dynkin index of the representation R, 1/2 for fundamental representation 
fermions, (N + 2)/2 for two-index symmetric representation fermions in the color group 
SU(N), and so on. Thus we expect to see multiples of 5 zero modes for sextet overlap 
fermions in our SU(3) case. 

However, ten years ago Heller, Edwards, and Narayanan discovered [58[] that the index the- 
orem applied to adjoint overlap fermions in background SU{2) gauge configurations failed: 
while 2T(R) = 4, they saw configurations with zero modes which were not multiples of 
four. More recently, similar results were reported by Garcia Perez, Gonzalez-Arroyo and 
Sastre[5^]. In simulations where the bare gauge coupling is large, I have seen configurations 
whose zero mode content was a not a multiple of 5. 

Fortunately, the authors of Ref. have studied the index theorem for sextet represen- 
tation fermions in quenched background SU(3) gauge field configurations and shown that 
in the continuum limit only configurations with multiples of five zero modes are found. The 
"fractional states" are apparently just a particular failure of the overlap action to capture 
topology when the gauge configuration is rough. 

This is not a problem for this project. The gauge configurations at the parameter value 
used to compute eigenvalues {(3 = 5.2) were smooth enough that all of the lattices I collected 
had Q = 0. I performed some trial simulations at a lower (3 of 4.8 and also saw only Q = 0. 
However, for these rougher configurations the cost of the overlap operator went up by about 
a factor of four, and it did not seem like a good use of computer time to continue running 
there. 

A map of the simulation region is shown in Fig. [TJ It is qualitatively similar to what 



we found with another bare action [llj. The line is the location in bare parameter space 
where the AWI quark mass vanishes, k c . The crosses show the location of the N t = 6 
deconfinement phase transition. To the left of this line, the system is confined and seems 
to be chirally broken. However, this region is little explored. It is unknown how (or if) the 
deconfinement line attaches to the k c line. To the right of the line, the system is deconfined 
on all observed volumes and spectroscopy shows parity doubling. This plus the smallness 
of the pseudoscalar decay constant leads me to conclude that the system is in a chirally 
restored state throughout that phase. All data used in this study comes from the weak 
coupling phase. 

Table fl] lists all simulation points use in this study. Since I will report results in terms of 
the AWI quark mass, I include it in the table. When there are several volumes, the quoted 
mass is from the smallest volume. However, generally the volume dependence of the quark 
mass is small ( a few digits in the least significant figure). In the analysis, data from any 
particular volume is plotted and used at its measured quark mass in that volume. As (3 
falls, the available parameter space in the deconfined phase shrinks to a smaller and smaller 
range of quark masses. 
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FIG. 1: Map of the bare coupling constant plane relevant to our sextet simulations. The solid 
line is the line of zero quark mass, k = k c . The crosses show the location of the confinement- 
deconfinement crossover at Nt = 6. 

III. FINITE-SIZE SCALING 

I begin this section with several pictures of spectroscopy to set the stage. Fig. [2] are plots 
of screening masses at /3 — 5.2 on 12 3 x 6 and 12 4 volumes. These two panels show that in 
the weak coupling phase the excitation spectrum is parity doubled and that the pseudoscalar 
decay constant becomes small as the quark mass vanishes. 

Throughout the weak coupling phase, masses depend strongly on the quark mass and 
weakly on the bare gauge coupling. Fig. [3] illustrates this feature for the correlation length, 
denned as the inverse of the pseudoscalar screening mass, showing 12 3 x 6 volumes on the left 
and 12 4 volumes on the right. Different plotting symbols correspond to different f3 values. 

Finally, I combine data from many volumes at one gauge coupling, (3 = 5.2, in Fig. HI 
The saturation of the correlation length at a scale proportional to the temporal length of 
the lattice is apparent. Superficially, this is just the Matsubara cutoff M ~ 2n/N t expected 
from the fermions' antiperiodic temporal boundary conditions. 

The relation between correlation length £ and quark mass given by Eq. [2] is only expected 
to hold when the system size L is much larger than £. When the correlation length measured 
in a system of size L (call it £5) becomes comparable to L, £l saturates at L even as m q 
vanishes. However, if the only large length scales in the problem are £ and L, then overall 
factors of length can only involve £ and L. For the correlation length itself, this argument 
says that 

£ L = LF(£/L) (14) 




TABLE I: Simulation points reported in this work. Since I am using the AWI quark mass as my 
independent variable, I catalog it along with the bare parameters (/?, k, volume). 



P 


K 


OiTflq 


volumes 


4.ZU 


n 1 ^7n 


n 1 1 n 

U.11U 


1Z X D 


A A n 
4.4U 


n 1 QQC 

U.looO 


U.lUb 


IZ X b 


a /in 
4.4U 


U.lo4o 


n n^i 
U.UOZ 


iz x b 


4.bU 


n 1 onn 

U.13UU 


n 1 r 

U.lo5 


12 x b, 12 


A c n 


n 1 onn 

U.132U 


U.Ubo 


12 x b, lz 


41. ou 


n 1 "a\ok 


n 097 
U.UZ I 


1 9 4 
IZ 


4.0U 


n 1 OQH 

U.lZoU 


n /i 1 
U.41Z 


1 o3 v > p. 1 o3 *✓ 1 o4 
IZ x b, lz x 0, lz 


A on 
4.8U 


n 1 Ten 
U.lzoU 


n onn 
U.oUU 


1 o3 c 1 o3 w q 1 o4 

lz x b, lz x 0, lz 


/i on 
4.8U 


n 1 o^in 

U.lzbU 


n oqo 
U.233 


1 o4 
lz 


/I on 
4.8U 


n 1 oTn 

U.127U 


n onn 

u.zuy 


1 o3 * / a 1 o3 * , 1 o4 

lz x b, lz x 8, 12 


A on 
4.8U 


U.1285 


f\ A Ad 

U.14b 


1 o3 w o. 1 o4 1 £?4 

lz x b, 12 , lb 


/i on 

4. 80 


n 1 onn 

U.129U 


n inn 

U.1U4 


1 o3 * - 1 o4 1 /^4 

12° x 8, 12 , lb 


A on 

4.8U 


n 1 onn 

U.13UU 


n nor 

U.U8b 


1 o3 v ^ 1 o3 1 o4 1 £?4 

12 x b, 12 x 8, 12 , lb 


5.20 


0.1100 


0.940 


12 3 x 6, 12 4 


5.20 


0.1150 


0.630 


12 3 x 6, 12 4 


5.20 


0.1220 


0.310 


12 3 x 6, 12 4 


5.20 


0.1235 


0.220 


12 4 


5.20 


0.1250 


0.195 


12 3 x 6, 12 4 , 16 4 


5.20 


0.1270 


0.096 


12 4 


5.20 


0.1285 


0.066 


8 4 , 10 4 , 12 3 x 6, 12 3 x 8, 12 4 , 16 3 x 8, 16 4 


5.20 


0.1290 


0.047 


12 3 x 6, 12 4 , 16 4 



where F(x) is some unknown function of A somewhat more useful version of this 

relation can be written by using Eq. [2] to say 

a = Lf(L*»m q ). (15) 

Then one can plot ^l/L vs L ym m q for many L's, vary y m , and look for the appearance of a 
smooth curve. The data from different L's will march across the x axis at different rates. 

A good data set to use is the one of Fig. HI (3 = 5.2. These are screening masses, so I will 
take L = N t regardless of the value of N s . A scan of Eq. [15] is shown in Fig. [5l Already one 
can see that a choice of y m ~ 1.5 pulls all the data from Fig. H]onto a single curve. 

To turn this observation into a number with an uncertainty is a bit awkward. We are not 
doing a fit, because the function f(x) of Eq. [T5]is unknown. Instead I will do the analysis 
in two different ways: 

First [6l|] take the full range of values of £l/L and slice it into a set of Nb bins. In the 
jth bin, define 

to) = ^ E m ^ m ( 16 ) 

idbin 

and 

i£bin 
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FIG. 2: Screening masses at (3 = 5.2 for two simulation volumes. Left: 12 3 x 6, right, 12 4 . Symbols 
show mass vs AWI quark mass: diamonds - pseudoscalar; squares - vector; octagons - axial vector; 
burst - scalar; crosses - pseudoscalar decay constant fps- 
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FIG. 3: Correlation length (inverse pseudoscalar screening mass) versus inverse AWI quark mass, 
on 12 3 x 6 (left) and 12 4 volumes (right). Plotting symbols are bursts, (3 = 4.2; crosses, (3 = 4.4; 
squares, (3 = 4.6; octagons, (3 = 4.8; diamonds, (3 = 5.2. 



Then define a goodness-of-fit parameter as 

Minimize this function with respect to y m and fold the whole procedure into a jackknife over 
all the different mass-L combinations to get an uncertainty. The division by (xj) 2 mimics 
the eye's attempt to minimize the fractional spread of the points about their average. Notice 
that bins with only a single entry do not contribute to x 2 . Binning the data introduces the 
possibility of a dependency on the number of bins. I assign an error by jackknifing the data 
set. In practice, if there are too many bins, there are seldom bins with more than one point, 
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FIG. 4: Correlation length (inverse pseudoscalar screening mass) versus inverse AWI quark mass 
at (5 = 5.2. Plotting symbols are for different simulation volumes, diamonds, 12 3 x 6; octagons, 
12 3 x 8; squares, 16 3 x 8; crosses, 12 4 ; bursts, 16 4 . 



and the jackknife error becomes large. 

A second approach follows the method of Ref . 62] . The idea is to use each data set (each 
different L value) to estimate f(x) and to find the y m which pulls the other L sets onto 
it. This is done inclusively; all data sets take a turn at being the fiducial. The quantity to 
minimize is 

p(y m ) = ^EEE f^P 1 - m m ^)) (19) 

^ over \ J-'j 1 

The interpretation of this long formula is that data set p is used to estimate the scaling 
function f(x). This is done by interpolation, either by polynomials or rational functions, 
using the recorded values of £l/L. The label "over" indicates that the sum only includes 
data from set j whose x values, Lj m mjj, overlap the range of x's of set p. The overall factor 
of l/N over counts the total number of points and guards against recording a zero value of P 
if there are no overlap. (Bhattacharjee and Seno actually consider powers other than two 
inside the sum.) P is minimized by the optimal y m . This method has an advantage over the 
first one that one does not need to bin the data, and a disadvantage that the interpolation 
algorithm has to be robust. This can be a problem for extreme values of y m . The number 
N over varies as y m is tuned. This could potentially make P discontinuous. 

The two methods produce similar results. At (3 — 5.2 there are 23 separate data points 
(values of m q and L), with C,l/L ranging from about 0.08 to 0.16. Asking for 5 bins gives 
(on average) 4 useful bins and y m = 1.43(25). Asking for 6 bins gives 5 useful bins and 
y m = 1.44(24). Asking for 8 bins produces on average 5 useful bins and y m = 1.51(38). 
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FIG. 5: Plots of £l/L vs m q L ym at = 5.2 for four choices of y m : (a) y m = 1.0, (b) y m = 1.4, (c) 
y m = 1.8 (d) y m = 3.0. Plotting symbols are for different simulation volumes, diamonds, 12 3 x 6 
(L = 6); octagons, 12 3 x 8 (L = 8); squares, 16 3 x 8 (L = 8); crosses, 12 4 (L = 12); bursts, 16 4 
(L = 16). 



With the second method, y m = 1.53(13) from a single-elimination jackknife over the data 
set. 

Bhattacharjee and Seno advocate taling an error from an approximation to the second 
derivative of P, 

/ P(y m (l + V )) Y 1/2 
Ay m = r]y m 12 In ^— 1 (20) 

With 77 = 0.1 this produces y m = 1.54(11). 

At (3 = 4.8 the situation is similar, but noisier. I again have four L's, from 12 3 , 12 3 x 8, 12 4 
and 16 4 volumes (L = 6, 8, 12, 16 respectively) and a total of 21 mass and size combinations. 
The data is noisier than the (5 = 5.2 data set. Fig. [6]shows the data before and after collapse 
to a line. A five-bin single jackknife fit gives y m = 1.21(32) while the second method gives 
1.41(26). The approximate second derivative gives y m = 1.40(20). A conservative summary 
of values and uncertainties is y m = 1.5(2) at (3 = 5.2 and 1.4(2) at (3 = 4.8. 

In this analysis, I made a particular choice of a geometry and a method of defining a 
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correlation function (screening masses, sensitive to the length of the antiperiodic boundary 
condition in the temporal direction). I do not believe that it is better than alternative 
methods one could try (for example, taking L 3 x L t volumes and varying L) . People interested 
in finite size scaling tests should try other geometries. 

In the statistical mechanics finite size scaling literature it is common to perform finite 
size scaling analyses on susceptibilities. I have experimented with this, but probably not 
extensively enough. Mocking up susceptibilities by integrating over correlation functions 
(x ~ J2x(0( x )0(®))) nas n °t produced interesting peaks, because the correlators - and 
their integrals - are dominated by short distance effects. This deserves more study. 



IV. THE CONDENSATE DIRECTLY FROM THE OVERLAP OPERATOR 

In this section and the next one, I only have data at (3 = 5.2. I ran at one k value, 
k = 0.1285. This corresponds to a small AWI mass. As I remarked, it was quite expensive 
to evaluate the overlap operator at stronger coupling. 

The condensate is measured in the usual way, with a vector of Gaussian random numbers 
defined on every site of the lattice, r/i, through an average over the random vectors, 

S K) = = ^HvlD(m q )-\ (21) 

i i 

and as usual for the overlap operator, D{m q )~ 1 is the subtracted, shifted Dirac operator 
Dijriq)- 1 = (D(m q )- 1 -l/(2r ))/(l-m q /(2r )). I accelerated the inversion by deflating with 
the eight smallest eigenmodes of D^D. This is all quite standard [52 ] . Even a small data set 
(twenty lattices, twelve random vectors per lattice) produces a nice signal (Fig. [7^), which is 
almost completely linear in the quark mass. This demonstrates rather dramatically that the 
gauge field configurations at one point in the weak coupling phase do not allow chiral valence 
quarks to form a condensate. The linear behavior is just the UV-dominated (proportional 
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FIG. 7: (a) Condensate and (b) S(m, mo) from valence overlap fermions as a function of the valence 
overlap mass, from background configurations at /? = 5.2, k = 0.1285, on 12 4 volumes. 



to m q /a 2 ) term in the expected expansion. I tried to remove it using a variation on a trick 
from Ref. j63|: Compute S(m 9 ) for each random number, for many masses simultaneously 
by using a multi-mass sparse matrix inversion algorithm. Then pick a mass as a fiducial and 
compute the difference S(m q ,m )i = T i {m q ) i — (m q /m )T>(m )i random number by random 
number. Finally, form the average over random seeds. The high correlations in the values of 
the condensate for different masses are removed from the error budget in S(m q , mo). I took 
nine masses ranging from am q = 0.01 to 0.25 and used arriQ = 0.01. The result is shown in 
Fig.[7)o. 

Unfortunately, a fit of S(m q ,m ) to a power, m^, produces a noisy result that a ~ 3. 
The exponent is unstable against the set of masses included in the fit. This is most likely 
(and one could not argue that it is not) just the nonleading UV term (m^loga). Thus 
this attempt fails. Something else is required, which is insensitive to the UV part of the 
condensate. That follows in the next Section. 



V. FINITE SIZE SCALING OF DIRAC EIGENVALUES 

To test the scaling law of Eq. I generated ensembles of lattices at one set of bare 
parameters, (3 = 5.2 and k = 0.1285, and a number of simulation volumes, and computed the 
lowest eight eigenvalues of the massless overlap Dirac operator on them. The bare parameter 
set was chosen to have a light valence quark mass and to be within the weak coupling phase. 
All of the configurations collected at this parameter value have zero topological charge. To 
see scaling, it is necessary to preserve the geometry of the simulation volume (to compare, 
for example L 4 lattices at different L's). My primary data set is L 4 lattices with L = 8, 10, 
12, and 16. The resulting eigenvalue spectrum is shown in Fig. [HJ By eye, it seems to show 
power law behavior. 

Fig [9] illustrates the quality of the data. The left panel shows simulation time histories 
of the lowest four eigenvalues of the 16 4 data set. Each measurement is separated by five 
HMC trajectories. The right panel shows the error on the average computed by blocking 
Nf, successive measurements together. The lowest eigenvalue is clearly the noisiest, but the 
autocorrelation time does not seem to be too large: I bin two successive lattices (iV& = 2 or 
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FIG. 8: Average values of 8 lowest eigenvalues of the sextet-representation valence overlap operator, 
vs 1/L, where the lattice volume V is defined to be V = L 4 . The actual volumes, moving from left 
to right across the graph, are 16 4 , 12 4 , 10 4 and 8 4 . 

At = 10 HMC units) together before averaging. 

Now I wish to extract an exponent from Fig. [HJ For theoretical input, all we have is 
Eq. [HJ the finite-size scaling formula. One does not know a priori if it applies to all the 
eigenvalues, only to the lowest eigenvalues, or if there is some minimum volume for which 
it applies. I will just proceed empirically: I will look at fits to individual eigenvalues, then 
groups of them. I will fit all the data sets or drop the smallest volume and fit only the larger 
ones. 

I begin by fitting individual eigenvalues (lowest, second, and so on) to a power law, 
ln(Aj) = Ai — plnL. I choose to fit to all four volumes, or the largest three. The individual 
data points in each fit are uncorrelated, of course. Fits and chi-squareds are shown in Table 
HT1 Examples of fits are shown in Figs. [TOlandfTTl 

I also fit groups of eigenvalues, using ln(Aj) = Ai—plnL for i — 1 . . . N eigenvalues. Now 
the data are correlated. I look at the quality of fits from uncorrelated fits, and then repeat by 
taking bootstrap averages of the data. The behavior is quite similar to the fits to individual 
eigenvalues. Some examples (with 2a bootstrap errors) on the average, chi-squared from 
uncorrelated fits: 

• Fit the lowest 4 eigenvalues and biggest 3 volumes: p = 1.59(5), x 2 jdof = 40/ (12 — 5) 

• Fit the lowest 4 eigenvalues 4 volumes: p = 1.54(4), x 2 /dof = 61/(16 — 5) 
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octagons for i = 2, diamonds for i = 3, and crosses for i = 4. 
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FIG. 10: Three-volume fits to individual eigenvalues, (a) the lowest eigenvalue (b) the first excited 
state (c) the second excited state (d) the third excited state. The volumes are (from the left) 16 4 , 
12 4 , and 10 4 . 
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FIG. 11: Four- volume fits to individual eigenvalues, (a) the lowest eigenvalue (b) the first excited 
state (c) the second excited state (d) the third excited state. The volumes are (from the left) 16 4 , 
12 4 , 10 4 and 8 4 . 

• Fit all 8 eigenvalues and biggest 3 volumes: p = 1.51(4), x 2 /dof = 126/ (24 — 9) 

• Fit all 8 eigenvalues and 4 volumes: p = 1.47(3), x 2 /dof = 180/(32 - 9) 



Examples of these fits are shown in Fig. 

The numerical value of the eigenvalues depends on the geometry of the lattice, but the 
scaling exponent does not seem to do so. I have data from 12 3 x 6 and 16 3 x 8 lattices. 
Fitting them as was done for the other data sets produces a similar exponent (similar results 
for individual or multiple eigenmodes, from a fit to the lowest four modes, p = 1.44(5)). See 
Fig. [13] for an example. With two volumes, of course, one can assume any scaling law that 
one wants. Nevertheless forcing a power law yields an exponent which is similar to that 
from L 4 volumes. 

At the end of this analysis I have many numbers, all rather similar but none really 
identical. What is missing from this section is some way of estimating the correction to 
scaling due to cutoff effects. In RMT, the larg er eigenmodes are the ones which are most 



affected by non-infrared physics. Kovacs [64j has pointed out that in distributions like 



p ~ \ a , the lowest eigenvalue has the broadest distribution and is most susceptible to finite 
statistics. Theoretical guidance is needed. Nevertheless, the observed exponents seem to be 
in quite good accord with the value of y m observed in Sec. 1111} y m ~ 1.5 at (3 = 5.2. 
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FIG. 12: Combined fits to several eigenvalues, (a) Three-volume fits to the 8 lowest eigenvalues, 
(b) four volume fits to the lowest 8 eigenvalues, (c) Three-volume fits to the lowest 4 eigenvalues, 
(d) four volume fits to the lowest 4 eigenvalues. The volumes are (from the left) 16 4 , 12 4 10 4 and 
in (b) and (d) 8 4 . 



TABLE II: Exponent p from fits to individual eigenvalues, from the largest three volumes (16 4 , 
12 4 , 10 4 ), or or to all volumes (add 8 4 ). 





3 volumes 




4 volumes 




mode 


P 


x 2 


P 


x 2 
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1.64(4) 
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1.60(3) 
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1.64(4) 
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1.59(3) 
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1.61(3) 
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1.56(2) 


16.5 


4 


1.58(3) 


18.7 


1.52(2) 


26.2 
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1.53(2) 
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1.49(2) 


19.2 
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1.51(3) 


10.0 


1.47(2) 


17.0 
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1.48(2) 
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1.45(2) 
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1.45(2) 
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1.43(2) 
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FIG. 13: Average values of 4 lowest eigenvalues of the sextet-representation valence overlap oper- 
ator, vs 1/L, for 12 3 x 6 and 16 3 x 8 volumes. As before, L 4 = V the lattice volume. The fit is to 
a common exponent, p = 1.44(5). 

VI. DISCUSSION 

The lattice-regulated system of Nf = 2 flavors of sextet fermions coupled to SU(3) gauge 
fields has a weak coupling phase which is chirally restored. Within that phase, hadronic 
correlation lengths (masses) show a weak dependence on the value of the bare gauge coupling. 
At one observation point in this phase, I find that the valence quark condensate vanishes as 
its (valence) mass is taken to zero. 

Motivated by the realization that slow running is very similar to no running (recall Eq. Hj), 
I analyzed the correlation length in the weak coupling phase as if its massless limit were 
critical, with the size of the relevant perturbation given by the AWI quark mass. I observed 
the collapse of correlation length data from many lattice sizes onto a common scaling curve. 
I did this at two bare parameter values. The observed exponent was identical within rather 
large errors. Roughly the same exponent governs the scaling of the values of low lying 
valence overlap Dirac eigenvalues at one of the couplings. My data cannot show whether the 
theory truly has only one relevant operator (with coupling proportional to the quark mass) 
and the gauge coupling is irrelevant, or if the gauge coupling is also running very slowly. 
This remains an open question. 

My analysis assumes that the correlation length would diverge at zero AWI quark mass 
in the infinite volume limit. The exponent y m is not unity, the value expected for a free field 
theory. The system thus has interacting dynamics. 
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Recall that y m = 1 — 7^. The perturbative expectation is 7^ = — 6C2(R)g 2 / (I6n 2 ) . 
Inserting the Schrodinger functional coupling quoted in the Introduction for these bare 
coupling values, perturbation theory predicts 7^ = —0.30 and -0.41 at the two couplings. 
This is a bit smaller than what I measure. (Of course, one could argue about the choice of 
a coupling in a perturbative formula.) 

Ryttov and Sannino j65| have a supersymmetric QCD - inspired beta function for gauge 
theories with higher dimensional representations of fermions. The anomalous dimension is 
predicted to be 

11C 2 (G) - 4T(R)N f 13 
lW 2T(R)Nf 10 K ' 

for the case studied here. My result disagrees with their prediction. 

This research leaves some obvious open questions. First, I have nothing to say about the 
strong coupling end of the weak coupling phase. 

Next, it is unknown if this theory has an IRFP or merely runs very slowly. Once that 
is known, results from this paper can be used for physics analyses. Even with the present 
statistical significance of my result, there are some rather definite conclusions which can be 
drawn once the answer to this question is known. 

If this theory does not have an IRFP, y m (g) can be used to run the condensate down 
according to Eq. [3J y m seems to be a bit larger than perturbation theory would give. Of 
course, collecting more points will give a nonperturbative determination of y m as a function 
(if any) of g 2 . 

If this theory is in fact an IRFP theory, my result indicates that it might not be a very 
spectacular IRFP theory, and it may not be phenomenologically interestin g. T his is because 
the exponent y m is small. Unitarity bounds for conformal field theories 66|, 67] constrain the 



scaling dimension of the leading scalar operator (which I am inferring will be the condensate, 
or at least what the fermions mass couples to) to lie in the range 3 > d — 4 — y m > 1. My 
y m = 1.5 is safely in that range. But in "unparticle" extensions of the Standard Model, new 
physics (NP) at scale M influences Standard Model (SM) physics at scale A through terms 
in an effective Lagrangian of the form 

^d NP +dsM 

C = + d °smO np (23) 

where the dimension of the new physics operator is d^p and the O's are operators in the two 
sectors. In the literature, d^p is generally desired to be as small as possible to enhance its 
effective coupling (the prefactor of the operator product). Luty and Okui j68|, for example, 
discuss models with d in the range 1-2. There is also an extensive literature relating the 
lower end of the conformal window to large y m , with = 2 perhaps having a connection 
with physics which closes the window. (See Refs. 69, 70, 71, 72, 73, 74|.) For sextet QCD, 



Nf = 3 is almost certainly close to the top of the conformal window, and Nf = 1 is almost 
certainly confined, so if Nf = 2 is conformal with y m ~ 1.5 there is not much of a story to 
be told. 

There are two other cases where y m has been measured. The first is SU (3) gauge theory 
with Nf = 16 fundamental representation flavors. Here, close to the top of the conformal 
window for its representation class, Hasenfratz [75| has found y m ~ 1. It is not surprising 
that y m is greater than 1 for Nf = 2 sextet QCD since y m is expected to increase with 
increasing distance from the top of the conformal window. The authors of Ref. 22| report 
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a determination of 7 m for SU (2) gauge fields and Nf = 2 adjoint fermions. They also find 
it is small. 

Understanding strongly-coupled beyond-Standard Model physics requires studying many 
theories and computing the anomalous dimensions for enough of them to be able to under- 
stand systematic trends. Finite size scaling studies are a useful way to do this. 
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